#include "aga.h"

#define LB_IMPLEMENTED 6
#define LBp  0
#define LBd  1
#define LBM2 2
#define LBM5 3
#define LBM6 4
#define LBHAC 5

#include <iostream>		// for std::cout
// #include <utility>                   // for std::pair
// #include <algorithm>                 // for std::for_each
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/connected_components.hpp>
#include <boost/graph/graphviz.hpp>

using namespace boost;

// create a typedef for the Graph type
typedef adjacency_list < setS, vecS, undirectedS,
        disallow_parallel_edge_tag > Graph;
typedef graph_traits < Graph >::vertex_descriptor Vertex;
typedef graph_traits < Graph >::edge_descriptor Edge;
typedef graph_traits < Graph >::out_edge_iterator EdgeIter;
typedef graph_traits < Graph >::degree_size_type degree_size_type;

Graph & generate_graph(GeneticData * data);

//! Massima cardinalita' di un pattern
int matrix_LBp(GeneticData * data);
int graph_LBd(std::vector < degree_size_type > &deg);
int graph_LBM2(std::vector < degree_size_type > &deg);
int graph_LBM5(std::vector < degree_size_type > &deg);
int graph_LBM6(std::vector < degree_size_type > &deg);
int hac(Graph g, unsigned int current_lb);

// void remove_nulls(GeneticData *data);
void remove_dominated_patterns(GeneticData * data);
void remove_singletons(GeneticData * data);
int is_included(int, int, GeneticData * data);
int remove_zero_row(GeneticData * data);
void delete_row(GeneticData * data, int r, int after, preprocT op_type);
int remove_zero_col(GeneticData * data);
void delete_col(GeneticData * data, int c);
// DEBUG
// void print_matrix (GeneticData *data, const char *str);

/* Preprocessing per riduzione delle dimensioni delle istanze:
 *
 * Le operazioni possibili sono:
 *  + eliminazione riga di zeri
 *  + eliminazione di colonna di zeri
 *  + eliminazione di un pattern <= in un altro
 *  + eliminazione di colonna di un singleton e della riga
 *
 * Non vengono considerate le sottoistanze dovute a non connessione.
 */

/*int lb_2(GeneticData *);
int lb_3(GeneticData *);
int lb_4(GeneticData *);
int lb_5(GeneticData *);
void bubble_sort(int *, int);*/

void preprocessing_initializer(GeneticData * data)
{
	data->stat.elaps[PREPROCESSING].start = clock();
}

void preprocessing_finalizer(GeneticData * data)
{
	data->stat.elaps[PREPROCESSING].end = clock();
}

//! Richiama le varie fasi di preprocessing.
/**
 * \param data	Puntatore alla struttura contente i dati dell'algoritmo
		genetico
 */
void preprocess(GeneticData * data)
{
	preprocessing_initializer(data);

// print_matrix (data, "Init");
	remove_zero_row(data);
// print_matrix (data, "Zero rows");
	remove_zero_col(data);
// print_matrix (data, "Zero cols");
	remove_dominated_patterns(data);
// print_matrix (data, "Dominated");
	data->connectedComponents = 0;
	remove_singletons(data);
// print_matrix (data, "Singleton");

	if (data->NumGenes > 2)
		find_mos_bound(data);
	else {
		data->MOS_Lower_Bound = matrix_LBp(data);
		data->connectedComponents = 1;
	}

	preprocessing_finalizer(data);
}

void find_mos_bound(GeneticData * data)
{
	Graph & g = generate_graph(data);
	std::vector < degree_size_type > deg(data->NumPieces);
	std::vector < degree_size_type > component(num_vertices(g));
	data->connectedComponents += connected_components(g, &component[0]);

	// genera ogni sottografo e calcola i bound

	for (Vertex v = 0; v < num_vertices(g); v++)
		deg[v] = degree(v, g);
	std::sort(deg.begin(), deg.end());

	int lbs[LB_IMPLEMENTED];

	lbs[LBp] = matrix_LBp(data);
//printf("lbs[LBp] = %d\n", lbs[LBp]);
	lbs[LBd] = graph_LBd(deg);
//printf("lbs[LBd] = %d\n", lbs[LBd]);
	lbs[LBM2] = graph_LBM2(deg);
//printf("lbs[LBM2] = %d\n", lbs[LBM2]);
	lbs[LBM5] = graph_LBM5(deg);
//printf("lbs[LBM5] = %d\n", lbs[LBM5]);
	lbs[LBM6] = graph_LBM6(deg);
//printf("lbs[LBM6] = %d\n", lbs[LBM6]);

	data->MOS_Lower_Bound = 0;
	for (int lb = 0; lb < LB_IMPLEMENTED - 1; lb++)
		if (data->MOS_Lower_Bound < lbs[lb])
			data->MOS_Lower_Bound = lbs[lb];
	lbs[LBHAC] = hac(g, data->MOS_Lower_Bound);
//printf("lbs[LBHAC] = %d\n", lbs[LBHAC]);
	if (data->MOS_Lower_Bound < lbs[LBHAC])
		data->MOS_Lower_Bound = lbs[LBHAC];
	delete & g;
}

/* Massima cardinalita' di un pattern */
int matrix_LBp(GeneticData * data)
{
	int m = data->NumPieces;
	int n = data->NumGenes;
	DataMatrix mat = data->DM;
	int max = 0;
	int pieces = 0;

	for (int i = 0; i < n; i++) {
		pieces = 0;
		for (int j = 0; j < m; j++)
			if (mat[i][j])
				pieces++;
		if (pieces > max)
			max = pieces;
	}

	return max;
}

/* min grado di un nodo + 1. Teorema 1 di Alice Miller */
int graph_LBd(std::vector < degree_size_type > &deg)
{
	std::vector < degree_size_type >::const_iterator i;
	for (i = deg.begin(); i != deg.end() && *i == 0; i++) ;
	if (i != deg.end())
		return *i + 1;
	else
		return 1;
}

/* Corollario 1 di Alice Miller: max(d[i] + 1 - 1), i=0..n */
int graph_LBM2(std::vector < degree_size_type > &deg)
{
	int max;		// deve essere intero con segno
	int first = 0;

	for (std::vector < degree_size_type >::const_iterator i = deg.begin();
	     i != deg.end() && *i == 0; i++)
		if (i != deg.end())
			first++;

	max = deg[first];
	for (unsigned i = 1; i < deg.size() - first; i++) {
		int lb = deg[i + first] + 1 - i;
		if (lb > max)
			max = lb;
	}

	return max;
}

/* Teorema 2 di Alice Miller */
int graph_LBM5(std::vector < degree_size_type > &deg)
{
	std::vector < degree_size_type >::const_iterator i = deg.begin();
	for (; i != deg.end() && *i == 0; i++) ;
	if (i != deg.end() && (i + 1) != deg.end()) {
		if (*i != *(i + 1))
			return *i + 2;
	}
	return *i + 1;
}

/* Corollario 2 di Alice Miller */
int graph_LBM6(std::vector < degree_size_type > &deg)
{
	std::vector < degree_size_type >::const_iterator i = deg.begin();
	for (; i != deg.end() && *i == 0; i++) ;
	if (i != deg.end() && (i + 1) != deg.end()) {
		if (*(i + 1) > *i + 1)
			return *i + 3;
	}
	return *i + 1;
}

/* void bubble_sort(int *array, int len)
{
	int swap = 1;

	for (int f = 0; f < len && swap; f++) {
		swap = 0;
		for (int l = len - 2; l >= f; l--) {
			int tmp;

			if (array[l] > array[l + 1]) {
				swap++;
				tmp = array[l];
				array[l] = array[l + 1];
				array[l + 1] = tmp;
			}
		}
	}
} */

int remove_zero_row(GeneticData * data)
{
	int i = 0;
	int &m = data->NumGenes;
	int &n = data->NumPieces;
	DataMatrix & mat = data->DM;

	for (int r = 0; r < m; r++) {
		bool all_zero = true;
		for (int c = 0; c < n; c++) {
			if (mat[r][c]) {
				all_zero = false;
				break;
			}
		}
		if (all_zero) {
			delete_row(data, r, -1, zero_row);
			r--;
			i++;
		}
	}
	return i;
}

void delete_row(GeneticData * data, int r, int after, preprocT op_type)
{
	int &m = data->NumGenes;
	int &n = data->NumPieces;
	DataMatrix & mat = data->DM;
	op_stackT & stack = data->preprocess_stack;

	if (stack.size == 0) {
		stack.size++;
		stack.ops =
		    (preproc_opT *) calloc(stack.size, sizeof(preproc_opT));
	} else if (stack.size <= stack.tos) {
		stack.ops = (preproc_opT *) realloc(stack.ops,
		                                    ((stack.size * 2) *
		                                     sizeof(preproc_opT)));
		stack.size *= 2;
	}

	stack.ops[stack.tos].op_type = op_type;
	stack.ops[stack.tos].original_index = r;
	stack.ops[stack.tos].operation.row_op.len = n;
	stack.ops[stack.tos].operation.row_op.original_row = (int *)(mat[r]);
	stack.ops[stack.tos].operation.row_op.insert_after = after;

	for (int i = r; i < m - 1; i++)
		mat[i] = mat[i + 1];
	mat[m - 1] =
	    (Pattern) (stack.ops[stack.tos].operation.row_op.original_row);
	m--;
	stack.tos++;
}

int remove_zero_col(GeneticData * data)
{
	int i = 0;
	int &m = data->NumGenes;
	int &n = data->NumPieces;
	DataMatrix & mat = data->DM;

	for (int c = 0; c < n; c++) {
		bool all_zero = true;
		for (int r = 0; r < m; r++) {
			if (mat[r][c]) {
				all_zero = false;
				break;
			}
		}
		if (all_zero) {
			delete_col(data, c);
			all_zero++;
			c--;
			i++;
		}
	}

	return i;
}

void delete_col(GeneticData * data, int c)
{
	int &m = data->NumGenes;
	int &n = data->NumPieces;
	DataMatrix & mat = data->DM;
	op_stackT & stack = data->preprocess_stack;

	if (stack.size == 0) {
		stack.size++;
		stack.ops =
		    (preproc_opT *) calloc(stack.size, sizeof(preproc_opT));
	} else if (stack.size <= stack.tos) {
		stack.ops = (preproc_opT *) realloc(stack.ops,
		                                    (stack.size * 2) *
		                                    sizeof(preproc_opT));
		stack.size *= 2;
	}

	stack.ops[stack.tos].op_type = zero_column;
	stack.ops[stack.tos].original_index = c;
	stack.ops[stack.tos].operation.col_op.len = m;
	stack.ops[stack.tos].operation.col_op.original_col = (int *)0;
	stack.ops[stack.tos].operation.col_op.insert_after = c - 1;

	for (int i = 0; i < m; i++) {
		int removed_element = mat[i][c];
		for (int j = c; j < n - 1; j++)
			mat[i][j] = mat[i][j + 1];
		mat[i][n - 1] = removed_element;
	}
	n--;
	stack.tos++;
}

//! Elimina i pattern inclusi in altri pattern.
/**
 * Elimina le righe che hanno uni inclusi in altre righe.
 * Ad esempio lo schema 0 0 1 1 e' dominato dallo schema 0 1 1 1 e pertanto
 * puo' essere eliminato.
 * In fase di postprocessing lo schema eliminato sara' da sequenziare dopo
 * lo schema che lo domina.
 *
 * \param data	Puntatore alla struttura contente i dati dell'algoritmo
		genetico
 */
void remove_dominated_patterns(GeneticData * data)
{
	int current, other;

	for (current = 0; current < data->NumGenes; current++) {
		for (other = 0; other < data->NumGenes; other++) {
			if (other == current)
				continue;
			if (is_included(other, current, data)) {
				// printf(" ** : %d <= %d\n", other, current);
				/* rimuove il pattern other */
				delete_row(data, other, current, subsumed_row);
				// printf(" ** Rimossa riga %d inclusa in %d\n", other, current);
				// print_matrix(data);
				other--;	// riparto dalla riga che ha sostituito other
				if (current > other)
					current--;	// aggiorno l'indice della riga di confronto
			}
		}
	}
}

/**
	Ritorna 1 se il primo schema passato e' incluso nel secondo.

	\param p1 L'indice del primo schema di taglio.
	\param p2 L'indice del secondo schema di taglio.
	\param data Contiene i dati sull'istanza.
	\return 1 se lo schema p1 e' incluso in p2, 0 altrimenti.
 */
int is_included(int p1, int p2, GeneticData * data)
{
	/* Suppone che p1 sia incluso in p2 */
	int j, included = 1;
//      printf(" ** Verifico %d <= %d\n", p1, p2);

	/* Per ogni pezzo in p1, verifica che sia anche in p2;
	   se cio' non avviene, allora p1 non e' incluso in p2. */
	for (j = 0; j < data->NumPieces; j++) {
		if (data->DM[p1][j] && (!data->DM[p2][j])) {
			included = 0;
			break;
		}
	}

	return included;
}

void remove_singletons(GeneticData * data)
{
	int &m = data->NumGenes;
	int &n = data->NumPieces;
	DataMatrix & mat = data->DM;
	int ii = -1;
	int jj = -1;

	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			ii = i;
			if (mat[i][j]) {
				// i,j possibile centro
				jj = j;
				// verifica che non ci siano "1" fino a fine riga
				while (++j < n) {
					if (mat[i][j])
						goto next_i;
				}
				// verifica che la colonna sia priva di "1"
				for (int k = 0; k < m; k++) {
					if (mat[k][jj] && k != ii)
						goto next_i;
				}
				delete_row(data, ii, -1, singleton);
				i--;
				delete_col(data, jj);
				j--;
				data->connectedComponents++;
				break;
			}
		}
next_i:
		;
	}
}

void postprocess(GeneticData * data)
{
	int &m = data->NumGenes;
	int &n = data->NumPieces;
	DataMatrix & mat = data->DM;
	op_stackT & stack = data->preprocess_stack;
	int **solution = (int **)(&(data->BestMosInd.Sequence));
	int k;
	bool stack_allocated = false;

	data->stat.elaps[POSTPROCESSING].start = clock();
	while (stack.tos) {
		stack_allocated = true;
		preproc_opT & op = stack.ops[--stack.tos];
		switch (op.op_type) {
		case zero_row:
		case singleton:
			m++;
			*solution = (int *)realloc(*solution, m * sizeof(int));
			for (int i = 0; i < m - 1; i++) {
				if ((*solution)[i] >= op.original_index)
					(*solution)[i]++;
			}
			(*solution)[m - 1] = op.original_index;
			for (int i = m - 1; i > op.original_index; i--)
				mat[i] = mat[i - 1];
			mat[op.original_index] =
			    (Pattern) op.operation.row_op.original_row;
			break;
		case subsumed_row:
			m++;
			*solution = (int *)realloc(*solution, m * sizeof(int));
			for (int i = 0; i < m - 1; i++) {
				if ((*solution)[i] >= op.original_index)
					(*solution)[i]++;
			}
			// cerca il punto in cui inserire il pattern rimosso
			for (k = 0;
			     k < m - 1
			     && ((*solution)[k] !=
			         op.operation.row_op.insert_after); k++) ;
			k++;
			for (int i = m - 1; i > k; i--)
				(*solution)[i] = (*solution)[i - 1];
			(*solution)[k] = op.original_index;
			for (int i = m - 1; i > op.original_index; i--)
				mat[i] = mat[i - 1];
			mat[op.original_index] =
			    (Pattern) op.operation.row_op.original_row;
			break;
		case zero_column:
			n++;
			for (int i = 0; i < m; i++) {
				int removed_element = mat[i][n - 1];
				for (int j = n - 2; j >= op.original_index; j--)
					mat[i][j + 1] = mat[i][j];
				mat[i][op.original_index] = removed_element;
			}
			break;
		default:
			exit(1);
		}
	}
	if (stack_allocated)
		free(stack.ops);
	data->stat.elaps[POSTPROCESSING].end = clock();
}

/*void print_matrix (GeneticData *data, const char *str) {
	int i, n;

	printf("%s\n%d %d\n", str, data->NumGenes, data->NumPieces);
	for (i = 0; i < data->NumGenes; i++) {
		printf("%3d: ", i);
		for (n = 0; n < data->NumPieces - 1; n++) {
			printf("%d ",
			       data->DM[i][n]);
			}
		printf("%d\n", data->DM[i][n]);
	}
	printf("\n");
}*/

//! Returns the max.
/*!
  \param a	A number
  \param b	Another number
  \return	The maximum
 */
static inline int max(int a, int b)
{
	if (a > b)
		return a;
	else
		return b;
}

//! Returns \f$\min\{a, b\}\f$.
/*!
  \param a	A number
  \param b	Another number
  \return	The minimum
 */
static inline int min(int a, int b)
{
	if (a < b)
		return a;
	else
		return b;
}

int hac(Graph g, unsigned int current_lb)
{
	degree_size_type deg, mindegree, graph_order, d1, d2;
	Vertex v1 = 0, v2;

	/*
	std::cout << "Initial graph" << std::endl;
	write_graphviz(std::cout, g);
	std::cout << std::endl;
	*/

	graph_order = num_vertices(g);

	mindegree = graph_order;
	for (Vertex v = 0; v < graph_order; v++) {
		deg = degree(v, g);
		if (deg == 0) {	// singleton da rimuovere
			mindegree = deg;
			clear_vertex(v, g);
			remove_vertex(v, g);
			graph_order--;
			v--;
		}
		if (graph_order == 0)
			return current_lb;
		if (deg > 0 && deg < mindegree)
			mindegree = deg;
	}

	current_lb = (current_lb > mindegree + 1) ? current_lb : mindegree + 1;

	while (true) {
		bool clique = true;
		for (Vertex v = 0; v < graph_order; v++) {
			deg = boost::degree(v, g);
			if (deg < graph_order - 1) {
				clique = false;
				break;
			}
		}
		if (clique) {
			current_lb =
			    (graph_order >
			     current_lb) ? graph_order : current_lb;
			return current_lb;
		}
		if (graph_order <= current_lb)
			return current_lb;

		// cerca i due nodi con grado minore e uniscili
		d1 = graph_order;
		v2 = graph_order;
		for (Vertex v = 0; v < graph_order; v++) {
			deg = boost::degree(v, g);
			if (deg < d1) {
				d1 = deg;
				v1 = v;
			}
		}
		if (d1 == 0) {	// singleton
			v2 = v1;
		} else {
			d2 = graph_order;
			for (Vertex v = 0; v < graph_order; v++) {
				deg = boost::degree(v, g);
				if (deg == d1) {
					for (Vertex w = 0; w < graph_order; w++) {
						Edge e;
						bool present;
						boost::tie(e, present) =
						    edge(w, v, g);
						if (w != v && present
						    && deg < d2) {
							v1 = v;
							v2 = w;
							d2 = deg;
						}
					}
				}
			}
			// merge v1 and v2 in v1
			if (v1 != v2) {
				EdgeIter edge_itr, end;
				for (tie(edge_itr, end) = out_edges(v2, g);
				     edge_itr != end; ++edge_itr) {
					Vertex n = target(*edge_itr, g);
					if (n != v1)
						add_edge(v1, n, g);
				}
			}
		}
		// remove v2
		clear_vertex(v2, g);
		remove_vertex(v2, g);
		graph_order--;

		/*
		write_graphviz(std::cout, g);
		std::cout << std::endl;
		*/
		mindegree = graph_order;
		for (Vertex v = 0; v < graph_order; v++) {
			deg = boost::degree(v, g);
			if (deg < mindegree) {
				mindegree = deg;
			}
		}
		if (mindegree >= 2)
			current_lb =
			    (current_lb >
			     mindegree + 1) ? current_lb : mindegree + 1;
	}

	return current_lb;
}

Graph & generate_graph(GeneticData * data)
{
	Graph *gp = new Graph(data->NumPieces);

	for (int pattern = 0; pattern < data->NumGenes; pattern++) {
		for (int piece = 0; piece < data->NumPieces - 1; piece++) {
			if (data->DM[pattern][piece] == 0)
				continue;
			for (int other_piece = piece + 1;
			     other_piece < data->NumPieces; other_piece++)
				if (data->DM[pattern][other_piece])
					add_edge(piece, other_piece, *gp);
		}
	}

	return *gp;
}
